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VARIANCE  REDUCTION  FOR  REGENERATIVE  SIMULATIONS,  Is 
INTERNAL  CONTROL  AND  STRATIFIED  SAMPLING  FOR  QUEUES 


by 


Donald  L.  Iglehart 
Stanford  University 
Stanford,  California 


Peter  A.  W.  Lewis 
Naval  Postgraduate  School 
Monterey,  California 


ABSTRACT 


We  discuss  two  methods  for  reducing  the  variance  of 
estimators  of  parameters  of  limiting  distributions  of  stable 
stochastic  processes  in  simulations.  The  methods  are  discussed 
in  the  context  of  the  simple  GI/G/1  queue.  Of  the  two  methods 
one,  which  we  call  an  internal  control  variable,  gives  a vari- 
ance reduction  which  is  roughly  uniform  over  values  of  the 
parameters  of  the  process  and,  in  particular,  works  well  for 
values  of  the  traffic  intensity,  close  to  1. 
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VARIANCE  REDUCTION  FOR  REGENERATIVE  SIMULATIONS,  Is 

INTERNAL  CONTROL  AND  STRATIFIED  SAMPLING  FOR  QUEUES 

by 

Donald  L.  Iglehart  and  Peter  A.  W.  Lewis 

1 . Introduction  and  Summary 

Simulators  are  frequently  faced  with  the  task  of  estimating 
a parameter  associated  with  the  limiting  distribution  of  a stochas- 
tic process  which  is  being  simulated.  A methodology  based  on 
regenerative  processes  for  obtaining  point  estimates  and  confidence 
intervals  for  such  parameters  has  recently  been  developed  in  Crane 
and  Iglehart  (1974a, b) , (1975a, b)  and  Iglehart  (1975),  (1976a, b) . 

In  this  paper  we  shall  indicate  two  techniques,  internal  control 
variables  and  internal  stratified  sampling,  which  might  be  used  in 
conjunction  with  the  regenerative  method  for  obtaining  additional 
variance  reduction  for  the  estimates.  To  illustrate  these  tech- 
niques we  shall  restrict  our  attention  in  this  paper  to  the  GI/G/1 
queue.  Other  applications  of  these  ideas  will  be  dealt  with  in 
future  publications,  as  well  as  uses  of  the  regenerative  property 
which  may  possibly  be  more  suited  for  obtaining  variance  reduction 
for  estimates  in  stable  stochastic  processes. 

In  the  GI/G/1  queue  we  assume  the  zeroth  customer  arrives 

at  time  tg  = 0,  finds  a free  server,  and  experiences  a service 

time  vQ.  The  n—  customer  arrives  at  time  t and  experiences 

a service  time  v . Let  the  interarrival  times  t - t = u , 

n n n-1  n 

nil.  Assume  the  two  sequences  {v  : n>0}  and  (u  : nil)  each 


consist  of  independent,  identically  distributed  (i.i.d.)  random 
variables  (r.v.'s)  and  are  themselves  independent.  Let  E{vn>= 
y 1,  E(un}  = X-1,  and  p=A/y,  where  0<A,y<°°.  Thus  y has 
the  interpretation  of  the  mean  service  rate  and  X has  the  inter- 
pretation of  the  mean  interarrival  rate.  The  parameter  p is 
called  the  traffic  intensity  and  is  the  natural  measure  of  conges- 
tion for  this  system.  We  shall  assume  that  p < 1,  a necessary 
and  sufficient  condition  for  the  system  to  be  stable. 

While  many  characteristics  of  interest  can  be  estimated 
using  the  regenerative  method,  we  shall  restrict  our  attention  to 
the  waiting  time  of  the  n—  customer,  W^  (time  from  arrival  to 
commencement  of  service) . To  obtain  a representation  for  the  pro- 
cess {W  : n>0}  let  = v -u  and  set  S =0,  S = X,  + +X  , 

n n n-1  n 0 n 1 n 

nil.  The  following  well-known  recursive  relationship  exists  for 

the  W ' s : 
n 

W0  ■ °'  Wn+1  ’ <Wl1+'  ni0- 


By  induction,  we  also  have 


W = max{S  -S,  : k = 0,l, ,n>,  n>0. 

n ri 


Using  the  strong  Markov  property  of  the  process  {Sn:  n>0) 

it  can  be  shown  that  there  exists  a sequence  of  integer- valued  r.v.'s 

(6^:  k>0)  such  that  8Q  = 0,  8^ < 8^+^,  and  wg  = 0 with  proba- 

k 

bility  one.  In  other  words,  the  customers  numbered  8^  are  those 
lucky  fellows  who  arrive  to  find  a free  server  and  experience  no  wait- 
ing in  the  queue.  The  fact  that  there  exists  an  infinite  number  of 
such  customers  in  the  GI/G/1  queue  is  a direct  consequence  of  the 
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assumption  that  p < 1 and  the  strong  law  of  large  numbers.  If 


we  let  “k = k>l,  then  represents  the  number  of 

customers  served  in  the  k—  busy  period  (b.p.)  and  they  are 

numbered  ^k-1'  ^k-l+1"  * " ^k-1^  * 

Next  define  the  sequence  {Y^:  ^-1)  by 

V1 

Y,  = I W . , 

k i=B  3 
3 pk-l 

the  sum  of  the  waiting  times  in  the  k—  b.p.  Since  the  queue 
is  stable  for  p < 1 we  have  Wn  =»  W as  n ■*■<*> , where  =»  denotes 
weak  convergence. 

Our  goal  is  to  estimate  E{w}  by  simulation. 

It  is  known  that  E{W}  < ® if  and  only  if  E{(x+)2}  <<=°.  We 

shall,  for  simplicity,  assume  that  E{v2}  <®  which  will  guarantee 

that  E{W}  <°°.  The  regenerative  simulation  method  is  based  on  the 

analytic  results  that  the  sequence  ^Yk,ak^  : k - ^ is  independent  and 

identically  distributed  and  E{W} = E{Y^}/E{a^} , the  ratio  of  two 

means.  This  suggests  using  the  ratio  of  estimates  of  E(Y.)  and 

1 n 1 

E(a..)  to  estimate  E (W)  . Thus,  if  we  let  Y(n)  =—  £ Y.  and 

_ 1 1 n n k=l  K 

a(n)  =—  £ a,  , where  n now  denotes  the  number  of  cycles  observed, 

n k=l  K 

then  a ratio  estimate  of  E (W) , for  example,  is 


W(n)  » Y(n)/a(n)  . 


(1.1) 


(In  the  sequel  we  drop  the  dependence  on  n unless  necessary.) 

Now  let  Zk  = Yk  ~ E^W^ak/  ^>1,  end  note  that  the  sequence 
{Z,  : k>l}  is  i.i.d.  and  E{Z.}  = 0.  We  assume  the  variance  of  Z.  , 


E{Z2}  » a2  = var{Yj{}  - 2 cov{Yk + var{aJ{}E2  {W}  (1.2) 

is  finite  and  positive.  Then  one  can  easily  show  that  as  n -*■ «® 

n1/2  [Y(n)/a(n)-E(W}]  ,, 

a/E{ap N(0,1),  (1.3) 

where  N(0,1)  is  a mean  zero,  variance  one  normal  random  variable. 
This  result  yields  a confidence  interval  for  E{w)  provided  we 
can  estimate  a/E(a^}.  A variety  of  estimates  have  been  studied 
and  are  reported  on  in  Iglehart  (1975a) . Here  we  just  mention  two. 
The  so-called  classical  estimate  for  a/E{a1>  is  given  by 

Si  = [Sn  - 2S12(Y/S)  +s22(yao2j1/2/5, 

where  Sn  is  the  sample  variance  of  the  Y^'s,  s22  of  the  ak's' 

and  Si2  the  sample  covariance  of  the  Y^'s  and  oi^'s. 

The  second  estimate  of  0/E{oii}  is  the  jackknife  estimate 
which  is  defined  to  be 

1/2 


S2  = ( I <e.-6)2/(n-l)  ]X/\ 
* i=l  1 


where  0 , = n (Y/a)  - (n-1)  ( T Y./T  a.),  l<i<n,  and  0 = i 7 0, 

j?*i  3 j^i  1 „ Ai=l 

is  called  the  jackknifed  estimator  of  E (W) . Both  Si  and  S2 

are  strongly  consistent. 

The  jackknifing  technique  is  known  to  work  particularly 
well  as  a confidence  interval  estimate  for  ratios;  for  a large 
number  of  cycles  n the  computational  problem  is  severe,  but  in 
that  case  the  technique  using  (1.3)  and  Si  works  well.  For 
details  see  Miller  (1974)  and  Iglehart  (1975)  . 


n 


w 


I 

i 

The  problem  addressed  in  this  paper  is  how  to  apply  variance 
reduction  techniques  with  the  ratio  estimator  W(n).  In  almost 
all  practical  situations,  where  in  particular  one  might  want  to 
compare  the  mean  waiting  times  of  two  different  queueing  systems 
(Iglehart  (1976)),  there  is  a premium  on  achieving  the  minimum 
possible  variance  for  the  estimation  in  the  given  computing  time 
(number  of  cycles  allowed) . Variance  reduction  techniques  for 
simulations  are  discussed  in  Kleijnen  (1974)  and  Gaver  and  Thompson 
(1973) , but  there  are  difficulties  in  applying  these  to  ratio 
estimates  and  regenerative  systems.  In  particular  variance  reduc- 
tion via  the  usual  control  variable  techniques  is  difficult.  A 
variation  of  this  technique,  which  we  call  internal  control  variables 
and  which  is  generally  useful  for  ratio  estimates,  is  introduced 
and  shown  to  give  considerable  variance  reduction  for  the  point  esti- 
mate W(n).  Another  technique,  internal  stratified  sampling,  is  also 
explored.  It  is  a natural  technique  to  use  but  appears  to  be  diffi- 
cult to  use  with  ratios.  Moreover,  it  becomes  less  and  less 
effective  as  p-*-l,  while  the  internal  control  technique  holds  up 
well  for  p close  to  1.  In  fact  the  internal  control  variables 
described  here  for  an  M/M/1  queue  give  a variance  reduction  which 
is  fairly  uniform  for  all  values  of  y and  ratios  p<  1.  Better 
results  can  be  obtained  for  particular  cases  of  the  parameter  p. 

It  will  also  be  apparent  after  the  development  of  the  vari- 
ance reduction  techniques  in  the  next  sections  that  the  two  tech- 
niques for  confidence  interval  estimates  discussed  above  apply  to 
the  estimates  after  variance  reduction  has  been  applied. 
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2 . Internal  Control 

Two  main  methods  for  variance  reduction,  antithetic  variates 
and  control  variables  (Gaver  and  Thompson  (1973))  have  been  put 
forward  for  use  with  queues  in  the  non- regenerative  situation.  Of 
these,  antithetics  has  very  limited  utility  beyond  the  simple  M/M/1 
situation  in  which  it  is  patently  clear  how  to  generate  two  reali- 
zations of  samples  which  give  negatively  correlated  estimates. 

That  scheme  for  generating  negatively  correlated  estimates  does 
not  work  with  the  regenerative  method  because  the  regenerative 
blocks  in  the  original  realization  of  the  queue  and  the  antithetic 
realization  get  out  of  synchronization.  No  alternative  way  has 
been  found  to  utilize  antithetic  variates  with  the  regenerative 
technique  and  we  feel,  along  with  many  computer  scientists,  that 
the  technique  has  limited  validity  in  systems  simulation. 

The  technique  of  using  a control  variable  and,  in  particular, 
a regression-adjusted  control  variable  (Gaver  and  Thompson  (1973) 
p.  587)  is  much  more  broadly  applicable  in  systems  simulation, 
although  it  is  again  difficult  to  adapt  to  the  regenerative  situa- 
tion. Briefly,  say  we  are  simulating  an  M/G/l  queue  to  estimate 
E (W)  with  the  non- regenerative  technique  of  averaging  the  first 
m waiting  times  to  obtain  an  estimate  of  the  unknown  E (W) , 


i m 

~ I 


W , 


j=l 


(2.1) 


The  same  random  numbers  used  to  generate  the  m non-exponentially 
distributed  service  times  are  used  to  generate  m exponential 
service  times  for  simulating  an  M/M/1  queue  whose  input  stream 
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is  identical  to  that  of  the  simulated  M/G/l  queue.  One  would 
then  expect  that  if  the  G-distribution  is  not  too  different  from 


the  exponential  distribution,  successive  waiting  times,  say  Wj , 

in  the  M/M/1  queue  realization  would  be  close  to  (positively 

correlated  with)  the  corresponding  waiting  times  in  the 

M/G/l  queue.  Consequently,  the  average  of  the  W^'s,  say  w^, 

will  be  positively  correlated  with  W and  one  can  form  a new 

m 

estimate 


w » W + B (W ' -E  (W* ) ) . 
mm  mm 


(2.2) 


Now  E (W^)  is  close  to  E{W' } = p/y (1-p)  for  m large,  so 

E (W  ) ~ E (W  ) and  the  variance  of  the  new  estimate  will  be  a minimum 
m m 


I i 


B = -cov(W  ,W‘)/var(W’ ) j 
mm  m 


(2.3) 


in  fact 


wm  s.d.(W  ) 
m m 


n»  = 1/2 


Wm  s.d.(W) 
m m 


= (l-r2)*'*. 


(2.4) 


where 


~ cov(W  , W ) 

r = corr  (W  ,W  ) = m 

m'  m o~ 

m m 


(2.5) 


There  are  a number  of  important  points  to  be  noted  about 


this  technique: 


(i)  It  allows  one  to  use  known  analytical  results  (such  as 
the  expected  value  of  the  limiting  waiting  time  in  an 


M/M/1  queue)  to  reduce  the  variance  of  a simulation. 
Such  use  of  analytical  results  is  a basic  principle  of 
simulation. 





(ii)  It  can  be  extended  to  non-linear  controls  or  to  multiple 
control  by  more  than  one  control  variable. 

(iii)  The  great  art  in  the  technique  lies  in  finding  a control 
whose  mean  value  is  known  and  which  is  highly  correlated 
with  the  estimator  which  is  being  controlled.  Thus  (2.4) 
says  that  |r|  must  be  close  to  0.9  to  reduce  the  stand- 

to  one  half  of  o~  and  this 


ard  deviation  o~  ^ 

w w 

m m 

generally  is  equivalent  to  reducing  the  required  sample 


size  for  a given  precision  by  a quarter.  Controls  which 
are  that  highly  correlated  with  the  estimate  are  not 
easy  to  find. 

(iv)  It  is  in  general  too  much  to  ask  that  the  correlation  and 
variances  in  (2.3),  (2.4),  and  (2.5)  be  known  analyti- 
cally. They,  therefore,  must  be  estimated  from  the 
simulation  data  and  this  will  reduce  the  variance  reduc- 
tion which  is  attained. 


Now  using  this  method  to  control  the  regenerative  estimate 
given  in  the  introduction , 


W(n)  = Y(n)/o(n)  , 


where  n refers  to  a fixed  number  of  cycles  in  the  queue,  one  runs 
into  the  same  problem  of  synchronization  as  with  the  antithetic 
case,  namely  n cycles  in  the  M/G/l  queue  may  take  a very  different 
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| 


f: 


t 

> 

\ 


'■  i 


a 


number  of  waiting  times  to  achieve  than  are  needed  for  n cycles 
in  the  M/M/1  queue.  Moreover,  the  correlation  between  and 

will  be  weak  and  made  even  weaker  by  the  use  of  the  ratio 
estimate.  This  problem  becomes  rapidly  apparent  in  simulation 
studies  of  the  technique. 

The  following  technique  which  we  have  called  internal 
(within  block)  control  has  been  developed  to  overcome  this.  It 
is,  in  fact,  a special  case  of  the  very  broad  technique  called 
concomitant  variables  in  Gaver  and  Thompson  (1973)  , p.  588  and 
will  be  illustrated  only  for  point  estimation  of  E (W)  and  E(a). 

Its  extension  to  other  quantities  discussed  by  Crane  and  Iglehart 
(1974a)  is  immediate  in  principle,  although  the  control  quantities 
discussed  below  may  be  different. 

2 . 1 Internal  Control:  Basic  Ideas 

The  idea  of  an  internal  control  variable  is  simple.  In  the 

estimate  W(n) , the  averages  Y (n)  and  a(n)  contain  n random 

variables  Y^  and  which  are  each  functions  only  of  the 

th 

interarrival  and  service  times  occurring  in  the  k — cycle  (or  b.p.) 

and  are  independent  of  the  other  interarrival  times  and  service  times. 

Thus,  it  is  natural  to  use  some  function  of  these  2a,  random 

k 

variables  to  control  each  Y,  and  a,  . 

k k 

The  naive  application  of  this  idea  is  that  if  we  can  reduce 
the  variance  of  both  the  numerator  and  denominator  Y(n)  and  a(n) 
we  will  reduce  the  variance  of  W(n),  but  we  will  show  that  the 
situation  is  more  complex  than  this.  We  will  denote  a function  of 
the  random  variables  in  the  k—  cycle  by  C(k)  but  will  generally 
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drop  the  index.  In  general  we  also  use  CT(k)  to  denote  a control 

for  the  numerator  (t  ;>)  in  the  ratio  estimator  and  Cg(k)  to 

denote  a control  for  the  denominator  (bottom) . 

Typically,  C(k),  or  simply  C,  could  be  the  difference 

between  the  service  time  v0  and  the  time  to  arrival  of  the  next 

6k 

customer  from  the  arrival  of  the  8.  th  customer,  namely 

^k+l 

This  difference  has,  in  a GI/G/1  queue,  a known  mean 
M 1 - X-'*'  and  large  positive  values  of  this  function  C(k)  corre- 
spond to  large  values  of  and  , and  vice  versa.  We  return 

to  specific  control  variables  and  their  computation  in  the  next 
sub-section . 

Note  that  one  can  control  either  the  top  or  bottom  of  the 
ratio  estimator,  or  both;  to  fix  ideas  assume  we  control  the  top 
and  have,  in  general,  an  internally  controlled  estimate 

n I + 8T(CT-E(CT) ) } 

WCT(n)  = — — 1 , (2.6) 

a 

where  gT,  as  in  the  usual  control  estimation  technique,  is  fixed 
so  as  to  minimize  var(WCT(n)).  In  practice  it  is  usually  esti- 
mated from  the  simulation  data. 

Now  the  quantity  a2/ E2{a},  where  a2=E{Zkl  is  given  at 
(1.2),  is  just  the  leading  term  in  the  asymptotic  expansion  for  the 
variance  of  the  ratio  estimator  W(n),  and  carries  over  to  the 
more  complicated  situation  (2.6)  to  give  (asymptotically  as  n -*•  °°) 


2 


where  Y'=Y/E{Y},  a’  =a/E{a},  and  = CT/E{Y} . Fora 
derivation  of  this  result  see  Cramer  (1946),  p.  354,  eq.  (27.7.3). 

It  follows  from  (2.7),  as  before,  that  to  minimize  the 
asymptotic  variance  of  W^T(n)  we  must  take 


-0T 


cov (Y ' -a ' ,C^,)  cov(Y,Ct) 


varlCp 


var  (Cn 


E(Y) 

E (a)  var (C  j 


(2.8) 


But  most  importantly  we  notice  that  C,£  must  be  highly  correlated 
with  the  difference  Y1  -a'.  To  achieve  this  is  much  more  difficult 
than  finding  a quantity  which  is  highly  correlated  with  either  of 
Y or  a,  simply  because  Y and  a are  highly  correlated  and 
increase  together.  In  particular  if  ct=l,  which  has  a high 
probability  if  p,  the  traffic  intensity  is  small,  then  Y = 0. 

Note,  too,  that  E (W)  =E(Y)/E(a),  the  quantity  we  are  trying 
to  estimate  in  the  simulation,  appears  in  (2.8)  for  the  top  control 
regression  coefficient.  Also  similar  equations  to  (2.7)  and  (2.8) 
pertain  to  the  case  where  the  bottom  of  the  ratio  is  controlled 
(in  this  case  a) , and  simultaneous  equations  can  be  derived  for 
and  8D  if  both  top  and  bottom  are  controlled.  The  additional 

X D 

complexity  in  estimating  BT  and  Bg  does  not  seem  to  be  justi- 
fied by  simulation  results  (discussed  later)  which  also  show  that 
if  only  one  control  is  used  there  seems  little  to  choose  in  terms 
of  reduction  achieved  between  putting  it  on  the  top  of  the  bottom. 


In  both  cases  the  control  must  be  highly  correlated  with  the 


The  above  results  are  generally  applicable  for  any  regenerative 
process  in  which  ratio  estimates  are  used.  The  choice  of  C is, 
of  course,  the  art  in  the  design  of  a simulation  with  variance 
reduction  and  is  considered  for  the  GI/G/1  and  specifically  the 
M/M/1  queue  in  the  next  section.  In  all  cases  this  choice  is 
limited  by  one 1 s ability  to  compute  analytically  E(C) . 

2 . 2 Internal  Control;  Design  Considerations 

We  discuss  here  the  design  of  internal  controls  for  the 
M/M/1  queue,  the  ideas  being  applicable  to  the  GI/G/1  case  with 
the  proviso  that  the  computations  might  be  considerably  more  diffi- 
cult. Simple  computations  of  E(C)  are  given  here  and  more  diffi- 
cult ones  in  the  Appendix;  we  do  not  distinguish  between  bottom 
and  top  controls,  since  both  must  be  correlated  with  the  difference 
Y-a,  and  we  drop  consideration  of  cycle  number,  since  all  variables 
are  within  the  cycles  which  have  identical  structure. 

Again,  we  are  considering  estimation  of  E (W) , but  of  the 
many  possible  controls,  those  listed  below  would  probably  work  as 
well  with  other  functions  of  W,  e.g.  percentiles.  The  controls 
are  listed  roughly  in  order  of  complexity  of  computation  and  of 
supposed  correlation  with  Y-a.  Thi:»  can  usually  only  be  guessed 
at  and  generally  the  more  elaborate  controls  which  might  have 
greater  correlation  with  Y-a  are  more  difficult  computationally. 

Superscripts  on  C are  labels  to  differentiate  the  controls. 

We  have  discussed  above  the  difference  X^  = Vq-u^,  whose  moments 
are  simple  to  compute.  Then  we  have 

i 

i 

I 
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if  a = 1 

(2.9) 

if  a > 2 . 

(2.10) 

It  is  easy  to  compute  E(C'  ')  for  M/G/l  queues  and 
possible  for  the  GI/G/1  queue.  Thus,  for  the  M/M/1  case  we 
have,  using  the  Markov  property  of  the  exponential  distribution. 


P{Vq<u^}  = P { a= 


a=l)  = f 
J0 

■ r 


tl-Fu(y) ]dFv(y) 


r*Y<u0-v  Y, 


e ( Vie  ^)dy  = 

0 i+p 


(2.11) 


which  goes  from  1 to  1/2  as  p=  X/y  goes  from  0 to  1;  furthermore, 

Pt“22>  - i-ITp  - ife  • (2-l2> 


Now  given  that  xi  = v0  “ ui  greater  than  zero,  the 

excess  v^  - u^  is  distributed  as  an  exponential  random  variable 
with  parameter  u.  Therefore 

E<X+)  = E(C(1>>  = 0 * 5^  + ^ i . (2.13) 

The  variance  of  can  also  be  computed. 

One  would  generally  like  to  obtain  more  correlation  of  C 
and  Y-a  when  Y is  large,  and  one  feels  this  can  be  done  by 
bringing  in  the  second  waiting  time.  Thus  we  have  as  control 
candidates 


/ 


c<2> 

= 0, 

if 

a 

= 1 

+ 4- 

(2.14) 

= xx  + X2 , 

if 

a 

2 2; 

C ( 3) 

= (X^+X2)+  = w2; 

(2.15) 

c(4) 

0, 

if 

a 

= 1 

x* + (x^+x2)+. 

if 

a 

(2.16) 

> 2. 

The 

control 

(2) 

Cv  ' is  WQ  = 0 if 

a = 

1, 

it  is  wq  + wi  if 

and 

it  is 

WQ  + Wx  + X^  if  a > 3. 

It 

is 

an  attempt  to  capture 

the  effect  of  the  first  two  waiting  times  without  the  additional 
computational  complexities  involved  in  computing  the  expectations 


of  ClJ'  and  c'  . 

Simulation  results  show  that  and  C ^ give  very 

(2 ) 

little  more  control  than  C for  which,  in  the  M/M/1  case, 

we  have 


E(C(2)  ) = 0 + EU'j+X*,^) 

= lTp{E(Xl+X2  I Xt>0} 

= E (xt  I xt  > 0)} 


r{i+E(X+  I X+  > 0)} 

■{f  + E<x2»}  =!55f4I¥f) 


(2.17) 


using  (2.12)  and  the  fact  that  X2  is  independent  of  X*.  We  use 
the  notation  E(x,A}  for  E{xiA>,  where  X is  a random  variable. 


A an  event,  and  1A  the  indicator  function  of  A. 

Similar  computations  go  through  for  and  C ^ 


from 


these  we  will  need  later  the  following  illustrative  result  (M/M/1 
case) . 


1 


--ye- 


(2.18) 


P{a-2}  = P{X*>0,X^+v1<u2} 

= P{Xi>0,v0-u1+Vl<u2} 

= T+p{P(u2-Y}  = l+p  * (y  + X) 2 = (l+p) 3 ' 


where  Y,  the  sum  of  two  exponential (y)  random  variables,  has 
a Gamma (y, 2)  distribution. 

From  this  we  also  get  that 


P{a>3) 


1 P 

l+p  (l+p) J 


P2 (2+P) 
(l+p) 3 ' 


(2.19) 


Now,  none  of  these  four  controls  is  specifically  designed 
to  be  correlated  with  the  difference  Y - a.  As  a result  they  work 
well  for  y and  X such  that  Y takes  on  values  much  larger 
than  a. 

To  fix  this  one  might  try 

C(5)  = c(l)  -a,  i = 1,2,3  or  4,  (2.20) 

since  the  mean  E(C^)  is  easily  calculated  if  E(C^)  is 
known  for  i = 1,2,3,  or  4,  and  E(ot)  is  known.  But  E(a)  is 
l/(l-p)  for  any  M/G/l  queue  (Cohen,  1969);  approximations  for 
the  GI/M/1  case  are  discussed  in  the  Appendix. 

There  is  an  additional  problem  of  dimensional  stability 
involved  in  using  C^,  as  with  , C^,  C^3*,  C*4  , in  that 

E(C^)  depends  on  both  y and  p,  while  E(a)  depends  only  on 
p.  Thus  control  is  not  uniform  across  the  whole  range  of  param- 
eter values. 
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To  avoid  this  one  use 


C(6)  = C(i)  - cl/m.  (2.21) 

This,  with  i = 2,  was  found,  in  simulation  studies,  to  be 

the  most  successful  control  variable  in  that  it  obtained  a variance 

reduction  which  was  uniform  in  y and  p (or  y and  X)  and  its 

mean  value  is  fairly  simple  to  compute.  These  simulation  results 

are  discussed  in  the  next  subsection. 

Note  that  multiple  control  variables  using  any  of  the  above 

controls  can  be  used.  In  particular,  one  need  not  take  the  differ- 
(2) 

ence  of  say,  C and  cl/m,  hut  may  use  a multiple  control. 

However,  the  fact  that  two  regression  coefficients,  say  6^^ 

(2) 

and  must  be  estimated  from  the  simulation  data  makes  the 

possible  gain  in  variance  reduction  of  dubious  value. 

Note  also  that  since  all  the  control  comes  from  within 
cycles,  there  is  no  reason  that  the  confidence  interval  estimation 
iques  referred  to  in  the  introduction  would  not  go  through 
for  t^.e  variance-reduced  estimates.  This  has  been  verified  in 
simulation  runs . 

2 . 3 Internal  Control;  Simulation  Results 

It  is  not  possible  to  verify  analytically  what  variance 
reduction  will  be  obtained  via  the  several  internal  controls 
listed  in  the  previous  section,  or  to  get  an  idea  of  the  magnitude 
of  the  effect.  Even  for  something  as  simple  as  it  is  diffi- 
cult to  compute  analytically  the  correlation  between  and 

Y-a  for  the  M/M/1  queue,  and  this  is  what  is  required  in  the 
equation  (2.5)  to  find  the  variance  reduction. 
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Thus,  we  resorted  to  simulations  to  verify  the  amount  of 


variance  reduction  obtained  and  the  relative  effectiveness  of  the 
various  controls.  In  the  final  simulations  all  runs  were  performed 
on  IBM  System  360/67  computer  using  the  LLRANDOM  package  (Learmonth 
and  Lewis  (1973))  which  generates  random  numbers  according  to  the 
scheme  given  by  Lewis,  Goodman,  and  Miller  (1969)  and  exponentially 
distributed  random  numbers  using  the  Marsaglia  "rectangle-wedge- 
tail"  method.  Tests  of  the  random  number  generator  are  given  in 
Learmonth  and  Lewis  (1974) . 

Of  the  extensive  simulation  checks  performed,  we  give  here 
only  a summary  of  the  conclusions  and  one  detailed  tabulation  and 
one  short  tabulation  in  the  case  of  the  most  suitable  control. 

(1)  The  controls  and  do  much  better  gener- 
ally than  C^,  with  little  improvement  over  obtained 

by  use  of  and  C^.  We  say  generally  because  results 

vary  with  \ and  y and  their  ratio  p. 

(2)  Subtracting  the  number  of  customers  served  in  a busy 
period  generally  improves  the  variance  reduction.  By  making 
it  dimensionally  stable  as  in  (2.21)  with  i = 2 we  obtain  a 


Table  1 shows  results  obtained  by  simulating  an  M/M/1  queue 


( 


I 


r. 


h 


■i 

) 

\ 

\ 


i 


i 


i 


i 


out  to  n = 2000  cycles  and  replicating  the  simulation  250  times  to 
estimate  the  variance  of  the  estimates  W(n),  WCT(n),  WCB(n), 
where  we  drop  the  n for  convenience.  Here,  we  have  specifi- 
cally that 


wCB(n)  " ~a 


1 n 

k ^ Yk 

nk=l  K 


nJi  tV8B(C<6)-E{c(6>))1 


where 


c<6)  = c,2> 


(2.22) 


The  estimated  precision  (standard  deviations)  of  the  estimates  of 
E (W)  are  given  in  brackets  under  the  estimates. 

The  results  in  Table  1 are  for  p = 0.5  and  three  values  of 
p;  the  results  are  not  very  different  at  different  values  of  p.  The 
case  p = 0.99  is  given  in  Table  2.  The  second,  third  and  fourth 
columns  in  the  Tables  give  correlations  between  the  control  and 
Y - a etc.,  from  which  the  theoretical  variance  reduction  can  be 
computed.  They  are  very  close  to  the  values  given  in  the  next 
to  last  column,  from  which  we  deduce  that  estimating  0_  and  3D 
affects  the  variance  reduction  only  slightly.  There  is  negligible 
effect  of  different  values  of  p for  fixed  p=0.5  and  fixed 
P = 0.99. 

Note  that  for  the  results  for  p = 0.99  given  in  Table  2,  the 
variance  reduction  is  73%  (about  the  same  as  for  p = 0.5).  For  the 
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0 
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case  where  the  control  is  on  the  top,  i.e.  for  Y,  the  variance 
reduction  is  not  quite  as  good  for  control  of  ot.  Note  too  that 
the  estimated  values  appear  in  some  cases  to  be  at  least  three  or 
four  standard  deviations  from  the  true  mean.  This  is  because  the 
estimates  W,  W and  W can  be  seen  from  the  100  replications 

to  be  non-normal.  In  other  words,  for  high  p(0.99),  the  simula- 
tion needs  to  be  taken  out  further  than  2000  cycles. 
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3.  Internal  Stratified  Sampling 

Another  technique  for  variance  reduction  which  can  be 
potentially  used  with  the  regenerative  method  is  stratified  sampling. 
For  a brief  description  see  Kleijnen  (1974),  p.  110.  In  essence 
this  uses  analytical  information  in  the  following  way. 

If  we  can  stratify  or  partition  a random  variable  X by 
its  values  (or  those  of  a concomitant  variable)  into  K strata 
labeled  k = l,2,...,K,  we  can  write  the  mean  of  X and  the  sample 
mean  X as,  respectively. 


K 

U = E(X)  = l p E(X  I Xe  strata  k)  ; 
k=l  K 


X 


, n 

- I 

r»  L. 


» a Xi 


* i 

str  1 


n 


+ l 

str  K 


n 


X. 


l 


(3.1) 


(3.2) 

(3.3) 


where  n^  = number  of  X^ ' s observed  in  strata,  and  pk  is  the 
probability  of  being  in  strata. 

Now,  if  the  p^ ' s are  known  and  we  substitute  them  in  (3.3) 
for  njc/n»  we  9et  Xst,  a stratified  estimate.  It  will  be  biased, 
since  the  divisions  of  the  sums  in  the  populations  are  random;  if 
the  numbers  observed  in  each  population  are  controlled  and  taken 
to  be  np^,np2»  etc.,  we  have  what  is  called  a proportionally 
sampled  estimate  X^g  with 


var(X  ) 
ps 


R 

var(X)  - l Pk(uk-tJ)  2/n, 


(3.4) 
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where  y^  = E{X  | X e strata  k}.  Thus,  because  of  the  use  of  prior 
analytic  information  we  always  get  variance  reduction. 

The  variance  of  Xgt  is  not  analytically  tractable,  but 
early  studies  reported  in  Kleijnen  indicate  it  is  close  to  (3.4) 
if  n and  all  the  n^ ' s are  sufficiently  large.  Our  simulation 
studies  with  stratifying  a have  confirmed  this;  because  of  the  ease 
of  computation  of  P { a=l } , P {a=2 } , P { a> 3 ) , as  in  section  2.3,  it  is 
natural  to  stratify  a.  Considerable  variance  reduction  is  obtained, 
especially  for  smaller  values  of  p,  the  traffic  intensity. 

Unfortunately,  when  the  quantity  a is  stratified  in  the 
bottom  (denominator)  of  the  estimator  W(n)  very  little  overall 
variance  reduction  is  obtained  unless  p is  small.  We  do  not 
understand  this  lack  of  variance  reduction  but  apparently  the 
correlation  between  the  stratified  version  of  a(n)  and  Y (n) 
is  reduced.  Analytical  studies  of  this  effect  are  very  difficult. 

It  is  also  possible  to  use  a stratified  estimate  of  a 


in  W(n)  and  to  control  the  top,  Y (n) , with  say  C 


using 


Cv  . This  works  well  for  small  p,  but  increases  the  variance 
as  p approaches  1.  Again  it  is  difficult  to  understand  this 
effect. 
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. Conclusions  and  Summary 

We  have  been  able  to  obtain  a worthwhile  variance  reduction 
using  internal  control  variables,  for  the  regenerative  estimate  of 
the  limiting  value  of  the  mean  waiting  time  in  an  M/M/1  queue. 

This  reduction  is  obtained  uniformly  over  all  parameter  values. 

It  is  fairly  certain  that  the  technique  will  work  well  with  any 
GI/G/1  queue  or  other  regenerative  stochastic  processes  or  systems. 
Internal  stratified  sampling  schemes,  however,  did  not  work  nearly 
as  well. 

The  techniques  can  be  extended  to  other  stable  stochastic 
systems,  such  as  the  Markov  chains  considered  in  Crane  and  Iglehart 
(1974b) . In  that  case  the  computation  of  the  mean  values  of  the 
controls  is  simpler  because  of  the  structure  of  the  Markov  chain. 

The  main  problem  in  applying  the  internal  control  variance 
reduction  technique  seems  to  lie  in  the  fact  that  the  estimator 
proposed  by  Crane  and  Iglehart  (1974a)  involves  a ratio  of  two 
random  variables,  and  these  are  difficult  to  work  with  in  general. 

An  alternative  which  will  be  considered  later  is  to  use  the 
existence  of  regeneration  points  more  specifically  to  obtain  vari- 
ance reduction  with  the  classical  estimator  W given  at  (2.1). 

m 

One  advantage  which  the  regenerative  estimator  W(n)  has  over 
Wm  is  the  ease  of  obtaining  confidence  interval  estimates  or  esti- 
mates of  the  precision  of  W(n)  and  W(n).  This  is  not  a draw- 
back if  the  simulation  is  large  and  more  than  one  (say  ten  or 
twenty)  realizations  of  the  queue  are  obtained. 

To  fix  ideas  note  that  we  can  write  W as 

m 
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W = ~ l w. 

m m . L •* 


V + Y ' ! 

k N(m)+lJ 


where  N (m)  is  the  number  of  completed  busy  periods  in  the  queue 
in  [0,m],  Y^  as  before  is  the  sum  of  the  waiting  times  in  the 
k—  cycle,  and  YN(m)+l  t*ie  sum  wa^t;*-n9  times  in  the  last 

incomplete  cycle. 

Now  it  is  possible  to  apply  internal  controls  to  each  Y^ 

in  the  sum.  Problems  arise  in  estimating  the  coefficients  8 in 

the  control  because  they  involve  a random  sum  of  random  variables. 

But  it  is  much  easier  to  find  a control  C for  Y.  rather  than 

k 

the  difference  Y^  - , and  also  it  is  still  possible  to  apply 

external  controls  as  well  as  internal  controls. 


These  ideas  will  be  followed  up  in  a later  paper. 


APPENDIX 


To  implement  either  the  internal  control  or  stratified 
sampling  techniques  certain  theoretical  parameters  associated 
with  the  GI/G/1  queue  are  required.  In  this  appendix  we  shall 
indicate  the  values  of  these  parameters  in  so  far  as  they  can  be 
calculated.  These  values  are  either  well-known  or  easily  calcu- 
lated. For  a reference  to  the  known  formulas  see  Cohen  (1969)  . 

We  begin  with  E{a^},  the  expected  number  of  customers 
served  in  a busy  period.  For  the  general  GI/G/1  queue  recall 

that  we  let  = v . - u and  S = X.  + . . . + X , for  n > 1, 
n n-l  n n 1 n 

with  Sq  = 0.  Then  a^  = inf{n>0:  SnSO}.  The  general  expression 
for  E(a^}  is  given  by 


E{a  } = exp{  l n 1P[S  >0]}, 
n=l  n 


an  impossible  expression  to  evaluate  in  general.  Another  useful 


expression  for  E{a^}  is 


E(a1)  = 1/P(W=0}, 


where  W is  the  stationary  waiting  time.  In  the  special  case  of 


M/G/l,  however,  we  have 


E {a^ } = (1+p)  . 


Now  for  the  queue  GI/M/1  we  can  use  (A.l)  and  the  stationary 
distribution  of  the  embedded  Markov  chain  to  conclude  that 


Etc^}  = ( 1-6 ) -x , 


/ 


where  6 is  the  root  inside  the  unit  circle  of 


z - U{y (1-z) } = 0 

-su 

with  U(s)  = E{e  },  Re  s>0,  and  u1  is  an  exponential  (p) 

r.v.  It  is  easy  to  check  that  6 = p for  M/M/1  queues.  Daley  (1975) 

has  recently  proposed  the  approximation  to  6 given  by 

S = a^ ( 1-p) 2 + 2 ( 1-b”1) p + (2b_1-l) p2 , 

where  a1  = P(u;L=0},  E{u1}  = 1,  and  b = E{u*h  This  approximation 

gives  good  results  in  a number  of  examples  calculated  by  Daley  (1975) 
and  may  be  useful  for  the  purposes  of  this  paper. 

Next  we  turn  to  the  computation  of  P{a1=l)  and  P{a^=2). 

For  the  GI/G/1  case  we  have 

Pla^l)  = PiS^SO) 

and 

P{a1=2)  = P{S1>0,S2£0}, 

both  of  which  can  be  worked  out  with  a little  effort.  For  the  M/M/1 

queue  , 

P{ai=l)  = U+p) 

and 

P{a-^=2}  = p(l+p)  \ 

For  the  M/G/l  queue 

P{a1=l>  = V( A) 

-Avn 

where  V ( A)  =E(e  ),  and  for  the  GI/M/1  queue 
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P{a1=l)  = 1 - U (M)  , 

where  U(s)  is  given  above.  For  the  M/E^/l  queue  and  the  E^/M/l 
queue  the  value  P{a^=2}  can  be  calculated  with  some  effort.  As 
these  expressions  are  cumbersome  they  shall  be  omitted. 


; ' i* 


Next  we  turn  to  various  partial  expectations  which  are 
needed  for  internal  control 

E{Sl+S2'al"2}  = EtS^S  >0}  + EiSjfS  >0} 


E^ai,Ctl22^  = ECa^}  - P{a1=l)  . 


In  the  special  case  of  the  M/M/1  queue. 


* i 


E{S1,S1>0}  = p/{y (1+p) }, 


Els2'sl>0>  - 12  (ife)Z  + TTFFF'  1 u'1' 


E(a1,a1>2}  = 2p/{  (1-P)  (1+p) } 
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